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The non equilibrium Green function formalism is today the standard computational method for 
describing elastic transport in molecular devices. This can be extended to include inelastic scattering 
by the so called self-consistent Born approximation (SCBA), where the interaction of the electrons 
with the vibrations of the molecule is assumed to be weak and it is treated perturbatively. The 
validity of such an assumption and therefore of the SCBA is difficult to establish with certainty. 
In this work we explore the limitations of the SCBA by using a simple tight-binding model with 
the electron-phonon coupling strength a chosen as a free parameter. As model devices we consider 
Au mono-atomic chains and a H2 molecule sandwiched between Pt electrodes. In both cases our 
self-consistent calculations demonstrate a breakdown of the SCBA for large a and we identify a 
weak and strong coupling regime. For weak coupling our SCBA results compare closely with those 
obtained with exact scattering theory. However in the strong coupling regime large deviations are 
found. In particular we demonstrate that there is a critical coupling strength, characteristic of the 
materials system, beyond which multiple self-consistent solutions can be found depending on the 
initial conditions in the simulation. We attribute these features to the breakdown of the perturbative 
expansion leading to the SCBA. 



I. INTRODUCTION 

Central to the field of molecular electronics are phe- 
nomena involving the interaction between the electron 
current and the internal degrees of freedom of the 
molecule investigated. In scanning tunnel microscopy 
(STM)i*2i£ the molecular vibrational modes (phonons) 
have been exploited to desorb or to move a molecule on 
a surface, paving the way for phonon assisted surface 
chemistry. At the same time, STM inelastic tunnelling 
spectroscopy uses the fingerprints of vibrations in the 
I-V curve to probe the orientation and/or to identity 
molecules on surfaces^^ 7 -. Switching devices exploiting 
phonons have also been reported 8 -. 

Broadly speaking, in molecular devices phonons are 
important for two reasons. Firstly, they play a role 
in transport^^ by opening new conductance channels 
through which the itinerant electrons can propagate, 
and by suppressing the transmission of purely elastic 
channels^. More dramatically, for large electron-phonon 
coupling the charge carriers become quasi-particles con- 
sisting of coupled electrons and phonons^. Secondly, 
from a technological point of view, phonons limit the effi- 
ciency of molecular devices because of energy dissipation. 
This causes heating, power loss and instability. 

Transport experiments at the nano-scale are difficult 
to interpret since the atomically precise device geometry 
is rarely known. Therefore one usually relics on atom- 
istic simulation techniques in order to understand the 
results. For clastic transport, when electron-phonon in- 
teraction is not considered and the electron-electron in- 
teraction is treated at the mean field level, methods of 
note for predicting the current flowing through devices 
include the non-equilibrium Green function formalism 
(NEGF^ii^OLi! and scattering theory (ST)^^^. 
Some of these methods have been adapted to include 



electron-phonon interaction, notably an extension of 
scattering theory (EST) n i 21 i 22 i 23 and the self consistent 
Born approximation (SCBA)^^ within the NEGF for- 
malism. In addition time-dependent methods for describ- 
ing correlated electron-ion dynamics have been recently 
proposed 2 ^. 

The focus of this paper is the SCBA. This is attractive 
from a practical point of view since it has moderate com- 
putational requirements and it has been used extensively 
for calculating transport properties of a number of differ- 
ent material system s 24 ' 25 . However, it is a perturbative 
approach appropriate only for weak electron-phonon (e- 
p) coupling. As the e-p coupling strength increases the 
SCBA will eventually breakdown, however it is unclear 
whether such breakdown is cither sharp or smooth with 
the e-p coupling strength. Our work explores this ques- 
tion in detail. 

The paper is organized as follows. We begin by pre- 
senting the NEGF formalism for a two probe devic o 27 ' 28 , 
and by recalling the foundations of the SCBA. We then 
consider a ID tight-binding model where the e-p inter- 
action in the scattering region is described by the Su- 
Schrieffer-Heeger (SSH) 29 i 30 Hamiltonian. The parame- 
ters for the model Hamiltonian are chosen for mimicking 
two systems which have been studied experimentally: H2 
molcules sandwiched between Pt electrodes (H2-Pt)2»2i 
and Au monatomic chains^ comprising R atoms (RC's). 
The parameters for H2-Pt are the same as those used 
by Jean and Sanvito^i, who previously employed exact 
scattering theory (EST) to describe phononic effects. In 
contrast to the SCBA, EST is valid for both strong and 
weak coupling and therefore it is a good benchmark for 
the SCBA. Accordingly we compare the SCBA results 
directly with EST over a range of different e-p couplings 
to establish the limit of validity for the SCBA and to 
investigate its breakdown. 
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II. METHODOLOGY 



A. Non Equilibrium Green Function Formalism 



A two probe device consists of two crystalline elec- 
trodes attached on either side of a scattering region, 
which is in general a collection of atoms breaking the elec- 
trodes translational symmetry. The leads arc also charge 
reservoirs, so that the device may be viewed as two charge 
reservoirs bridged by the central region. Thcrmodynam- 
ically we characterise the left-hand side (L) and right- 
hand side (R) lead by defining their chemical potentials 
/iL&nd /xr. If /iL = /jr. equilibrium is established and no 
current flows. When fiL ^ Mr the system is dragged out 
of equilibrium, and net charge will move from the reser- 
voir with the higher chemical potential across the central 
region to the reservoir of lower chemical potential in an 
attempt to re-establish equilibrium. If a battery is at- 
tached to the two reservoirs keeping ml — MR = °V (V 
is the bias and e the electron charge) the system cannot 
return to equilibrium and will eventually reach a steady 
state with a constant current flow. 

At the Hamiltonian level the problem can be formu- 
lated by using a basis set comprising a linear combination 
of atomic orbitals (LCAO). It is convenient to write the 
Hamiltonian of the semi-infinite periodic leads in term of 
principal layers (PLs) 27 ' 32 . These are cells that repeat 
periodically and constructed in such a way that the in- 
teraction between PLs extend only to nearest neighbours 
(see figure [1]). Thus the N x N matrices H\ and Hq 
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FIG. 1: Schematic representation of a system composed of 
two semi infinite leads and a scattering region (rectangular 
dashed box). The matrices Ho and Hi describe the lead prin- 
cipal layers, Hm describes the scattering region and -Hlm, 
-Hrm the interaction between the scattering region and the 
last principal layers of the leads. 



describe respectively the interactions between PLs and 
within a PL. The scattering region in general is described 
by M basis functions. The Mx M matrix, Hm, describes 
its internal interaction, while the matrices -Hlm (N x M) 
and /Trm {M x N) contains the interaction between the 
PLs of the leads adjacent to the scattering region and 
the scattering region itself. The entire system is thus 



described by the infinite tri-diagonal Hamiltonian TL 
( 
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Time reversal symmetry sets H—% = H\, Hml = H L 



LM> 



and Hmr = H^ M . The retarded Green function, Q 11 
associated to the entire system (leads plus scattering re- 
gion) is defined as 



(i) 



where u>' = lim,5^o+ ( w + u IS the energy and T is 
the infinite dimensional identity matrix4^ 

For transport calculations however one does not need 
the Green's function of the entire system but only that 
relative to the scattering region, Gm, in presence of the 
leads. This can be written as^ 



G M (w) = [w'J] 



M 



H 



M 



E L (w) - E R (w)]" 



(2) 



where the presence of the leads have been accounted via 
the introduction of the self-energies for the left- and right- 
hand side lead El (a;) and Sr(w). Im is the M x M 
identity matrix. The self energies are M x M matrices 
defined as 



Sl — Hml 9l Hlm , Sr — Hmr 9r Hrm 



(3) 



where <7l(^) and <7r(w) are the retarded surface Green 
functions of the leads, namely the retarded Green func- 
tions of the isolated semi-infinite leads evaluated at the 
PLs adjacent to the scattering region. These are calcu- 
lated by considering the retarded Green function of the 
corresponding infinite system (periodic) and by applying 
appropriate boundary conditions^. External bias volt- 
age is introduced under the assumption that the leads are 
good metals maintaining local charge neutrality. The ef- 
fect of a bias is therefore only that of shifting rigidly in 
energy the leads electronic structure, so that 



SL/R(w,y) = S L/R (a;±eV/2,0) 



(4) 



We now proceed to evaluate the non-equilibrium 
charge density in the scattering region and the two-probe 
current by using the NEGF scheme^ 2 -. The lesser (<) and 
greater (>) Green functions G m (uj) are defined as 



with self-energies 



ck=L,R 

i [ng(w) - l]r a (w) . 



(5) 



(6) 

(7) 

(8) 
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Here np(u>) = np(uj — /i Q ) is the Fermi function evaluated 
at w - /i„ and temperature T, /i Q = i5p ± with 
E-p the leads Fermi energy, and we have introduced the 
coupling matrix for the a lead (a = L/R) 



r a (w,VO=i[Ea(w,V)-Ei(w,V)] 



(9) 



The non-equilibrium charge density matrix for the scat- 
tering region is 



1 

2tt7 



dw'G<( W ') 



(10) 



If Hm has a functional dependence on p the equations @ 
and dTUJ) can be solved self-consistently. The net current 
flowing through the device is then 

2 poo 

Jun P (V) = — J duj TrTLG^rRGMK^F - n*) , (11) 

where the subscript "unp" stands for "unperturbed", 
meaning that no e-p interaction is included. The term 
T(ci>, V) = TtPlGi^PrGm] is the standard Landaucr 
Biittikcr transmission coefficient, although in this case 
it is explicitly bias dependent. The conductance G in 
the linear response limit is 



G 



2£ 
h 



T(E F ,0), 



while more generally at a given bias V one has 



G(V) 



dJ u 



dV 



(12) 



(13) 



B. Self Consistent Born Approximation 

We now discuss the main concepts associated with 
introducing e-p scattering into the NEGF transport 
scheme. In general, inelastic scattering produces loss of 
phase coherence, similarly to what happens when an elec- 
tron is absorbed by a reservoir. In fact one may think 
of inelastic processes as resulting from the coupling of 
the scattering region to a "fictitious" charge reservoir— 
that does not exchange a net current. Thus c-p interac- 
tion can be introduced via a self-energy E p h(w) and the 
retarded Green function becomes 

GmH = [a;'7 M - J ff M -E L (a;)-S R (a;)-S ph (a;)]- 1 . (14) 

The exact form for E p h(w) is unknown, however conve- 
nient approximations can be derived from the perturba- 
tive expansion over the e-p coupling strengt h 16 i 24 i 25 i 28 i 34 . 
In this work we consider the SCBA where only the 
Hartree and Fock diagrams of the perturbative expansion 
are retained (see figure [2]). This is equivalent to evaluat- 
ing the first order diagrams at the interacting electronic 
Green function. Thus the phonon self-energy reads 



E ph H = E F (co) + v/ 



(15) 



where the retarded Hartree (H) and Fock (F) contribu- 
tions to the self-energies are respectivel y 28 ' 35 

SH = 'Erf ^ A Tr[G<K)A^],(16) 



e f H = £ 



[E>( W )-E<( W )] 



A L 

i 

~2 



1 

2^ 

?0{E>(u/)-E<(u/)}(a;) 



(17) 



In equation (fTT|) H^ 1 is the Hilbert transform 

n x {m}(y) = -v [ X dx- f{ 

7T J-oo : 



(18) 



x-y 

and V stands for the principal part of the integral. The 

M 

| Hartree ^Fock 
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FIG. 2: Diagrammatic representation of the Hartree- Fock ap- 
proximation, (a) the self-consistent proper self-energy (the 
shaded circle) is obtained from the first-order Hartree- Fock di- 
agrams evaluated using the interacting Green's function (dou- 
ble line). This is equivalent to re-summing all the diagrams^ 
in (b). 

phonon energy and e-p coupling matrix for a particular 
mode A are respectively il x and M x . Finally the e-p 
lesser and greater self-energy are given by 



(19) 



E>M = M> 



(N x + l)G>{Lu±n x ) 



N x G>{ut^x) 



M 



(20) 
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which is simply a sum of the lesser self-energies over all 
the possible modes A. The occupancy of each phonon 
mode is N\. 

We assume that the phonons are in thermal equilbrium 
with a bath so that for a given temperature N\ is simply 
the Bose-Einstein distribution N\ = (e !!j / ftT - 
with Kb the Boltzman constant. From equation (f2T))) we 
note that the lesser (greater) self-energy contains con- 
tributions from two scattering processes: electrons with 
energy lu — fl\ may absorb (emit) a phonon and/or elec- 
trons with energy to + fl\ may emit (absorb) a phonon 
of energy Q \. When T ~ 0, electrons may only emit 
phonons since N\ ~ 0, i.e. no phonons are present in 
the scattering region provided that the phonon lifetime 
is much smaller than that of the electrons. The total 
lesser self-energy of equation ([5]) must be adjusted to in- 
clude the phonon lesser self-energy, 



J2 E°*( W )+£>( W ) 

q= L,R 



(21) 



and G>(u>) from equation ([5]) are now evaluated using 
the perturbed Green function and lesser self-energy from 
equations (fT4|) and (|2~T|) . The general expression for the 
interacting (including c-p coupling) current^ through 
lead P may be written as the sum of an elastic and an 
inelastic contribution 



where 



j on = 4(v) + jUv), 

4 = T / dtJ Tr[r /3 G M r Q G'L](n^ - n§) 

" J— oo 



(22) 



and 



T0 _2e 



duj Tr 



(23) 



III. NUMERICAL METHOD 

A. Model Hamiltonian and Coupling Matrices M x 

The systems under investigation are ID linear atomic 
chains described by a s-orbital nearest-neighbour tight- 
binding model. The scattering region comprises R atoms 
plus one PL (one atom) from each lead, so that it contains 
M = R + 2 orbitals. Henceforth we refer to this system 
as RC. Furthermore, we assume that the two leads are 
identical. The matrices H , Hi for a ID tight-binding 
model reduce to c-numbers: el = Hq and 7l = Hi, where 
£l j 7l are the lead onsite energy and hopping parameter 
respectively The leads' Hamiltonians thus read 

-1 -2 

= e L C k+7L ]T [4a + i + 4 +1 a} , (24) 

i— — oo i— — oo 

oo oo 

Ur = e L c * Ci + 7L E Cl + l + c Ui c i\ '( 25 ) 
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FIG. 3: Schematic diagram of the simple monatomic systems 
considered here. It is composed of two semi infinite leads and 
a scattering region. The scattering region is marked with a 
dashed rectangle. Inelastic scattering is effective only in the 
scattering region. 



where {cj, Cj} are the electronic creation and annihilation 
operators at site i. The interaction Hamiltonian between 
the leads and the scattering region is 



Hlm = 7lm[cLiC + cjc_i] 



(26) 



Hrm = 7LM[Cfj +1 CR + 2 +c T R+2 c R+ i] , (27) 

where for our setup 7lm = 7l- 

For a ID system the lead self-energies are analytical 



where 



A(w) = r 



e£(h)M = A(w)-i 



,.rV) 



larl < 1 



x - Vx 2 - 1 [acj > 1 ' 



r(w) - -2r fl(i- |x|)(Vi-x 2 ) 



(28) 

(29) 
(30) 



and Tn = 



7lm 
7l 



- £ l) 
27l 



Finally e-p interaction is included into the scattering re- 
gion in the form of an SSH Hamiltonia n 29 i 30 , comprising 
three terms 



Hm — H c + H cp + H p h 



(31) 



with 



R+l 

H c = 5>-t 

i=0 i=/Lj 



^U + J2iU c l c J + c h)> ( 32 ) 

H ph = Y(b{b x + hhw x , (33) 



A max 



i=R+2 



i=R+2 



H *p = EE M ^ ( 6 a + &A)(ct Cj + cjc). (34) 

A=l 

is the electronic Hamiltonian of the scattering re- 
gion with onsite energies and the unperturbed hop- 
ping parameters 7?- (we assume 7° = for j ^ i ± 1). 
i/ph is the non-interacting phonons Hamiltonian, writ- 
ten in terms of the phononic creation and annihilation 
operators {b\,b\} and the phonon energies tt\ = huj\, 



with the index A running over all the modest. The final 
term, H cp , is the Hamiltonian describing the e-p inter- 
action within the scattering region. The details of such 
interaction are included in the e-p coupling matrices M A - . 

In order to calculate the matrices M A - and the longi- 
tudinal phonon frequencies we consider a simple nearest 
neighbours elastic mode l 23 ' 24 i 38 in which 



M;. 




(35) 



In equation (|35[) the orthonormal vectors e\ represent the 
ionic displacement associated to each mode A, m, is the 
mass of the atom at site i, and the constants oiij are the 
e-p coupling parameters. These latter are defined as the 
first order coefficient of the expansion of the tight-binding 
hopping parameter 7^ about the atomic equilibrium po- 
sitions 



Hi 



1% + OLij (t 



(36) 



where Ui is the displacement vector of the atom at site 
i. From equation (|36[) it follows that = —ctji. In 
the nearest neighbour approximation the interaction is 
restricted to electrons moving between sites i and i ± 1. 
Finally, we note that the eigenvectors e\ are real. This 
implies that the matrices My are real and symmetric 
with non-zero matrix elements for i = i ± 1, so that Mh 
has a tri-diagonal form for longitudinal phonons. We 
note that, although it is possible calculate the coupling 
parameters using first-principles electronic structure 
methods, here we set ay = a and a is taken as a free 
parameter. 



Initialisation 



Numerical Parameters 
Hartree:E B , cp 1; cp^, rp x , Poles 
(dE) H ,t H 

Fock:Np,(dE) F ,t F 
Bias: V init ,V fmal ,dV 
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M\ ^[kcke^m] 
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FIG. 4: Scheme of our numerical procedure for self- 
consistency. 



B. Numerical Integration and Self Consistency 



where the ratio matrices i? A and their weighting coeffi- 
cients r\ are given by 



The flowchart in figure 2] outlines the numerical pro- 
cedure used to calculate the interacting current J^(V) 
and the differential conductance G(V). Each simulation 
can be partitioned into three steps. The first two consist 
of two self-consistent loops which calculate the phonon 
self-energies E H and S F respectively. These are used in 
the third step to evaluate J^{V) by using the equations 
(J22J and (PI) . 

Let us now discuss the three steps in some detail. We 
start by writing E H as an explicit function of the density 
matrix p 



S H (P) 



-4^Tr 



'nl 



M 



(37) 



where we assume that all the elements of are inte- 
grable. S H is thus nothing but a weighted sum of the 
matrices M x . This can be written in the form 



M 

E 

A=l 



(38) 



R> 



'a 



17° 



|M A | 



mm 



(39) 



M-l 



4|M A | max x - 



fi A | 7 °| 



1=1 



For a given mode A, the largest matrix element of the 
matrix M A is denoted as |M A | max . |7°| m i n is the smallest 
among the hopping parameters of the unperturbed sys- 
tem (no e-p coupling), which, by construction, is equal 
to |-R A |max- The matrices R x are independent of the e-p 
coupling a and simply reflect the symmetry of the specific 
phonon mode considered. Thus r\ measures the maxi- 
mum fractional modification of the elements 7?. in the 
electronic Hamiltonian as the result of e-p coupling. 

The first self-consistent loop of figure [4] begins by 
choosing initial values for the weighting coefficients, 
{fAj-init = { r i i r 2i which are used to construct 

the density matrix at the first iteration, p 1 . Then both 
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{r\} and p are varied until the convergence condition 



„n-l| 



„n-l I 



(41) 



is met for the chosen tolerance in- The density matrix 
p is obtained by integrating G < as in equation (jTOJ) . In 
order to perform this integral we first take E F = under 
the assumption that this has little effect of the conver- 
gence of p (first self-consistent loop). Following previous 
work s 27 ' 39 we write p = p cq + py , where 



duj ImfGjvi] n F 



and 



1 f°° 

Pv = 7T duj GmWIa {nf - 4) 

^ J-oo 

At equilibrium (p^ = /iR = fx) one has p = p cq . 
a) 



(42) 



(43) 




to iterate {E F ,E^ h } and {Gm,G^} until a second con- 
vergence condition is met 



^Maxj|[Gft(w<) 



7V F 



G^T 1 ^)]! < *f.(44) 



Note that the condition is over the largest of the ma- 
trix elements and runs over the iV F energy points u>i of 
the entire grid. Note also that the tolerance used, tp, 
is in general different to that used for the Hartree term. 
The Hilbert transform required for calculating the imag- 
inary part of E F (eq. (jTTJ)) is done by using a convo- 
lution method combined with a fast Fourier transform 
algorithm^. In order to avoid end-point corrections we 
choose a grid of sufficient range while the grid fineness 
(dE)p must be sufficiently fine to resolve phononic fea- 
tures which lie in the meV range. 

Table [J shows the numerical and system parameters 
used in our simulations. The parameters for the H/2-Pt 
junctions are identical to those used by Jeanii within the 
EST treatment of phonons. This set produces the same 
unperturbed G(0) ~ .97Go- The spring constants and 
masses are chosen to give longitudinal phonon modes of 
energies 63 meV (CM mode) and 432 meV respectively 
while the ratio matrices for these two modes are 



E F -V/2 E F +V/2 



3.2 























-3.2 





-3.2 






(45) 



FIG. 5: Schematic representation of the integrals in equations 
(|43l) and (|42l) respectively. In a) the integration of pv is bound 
between Ef + eV/2 and Ef — eV/2 by the Fermi functions, 
rpi is the number of points of the real axis energy grid. For 
p cq in b) the number of grid points along the circular path 
(epi) and the path in the complex plane parallel to the real 
axis (cp2) must be chosen, as well as the poles in the Fermi 
functions. 

As shown in figure [5^, the integral of py (Eq. (|4"3"| ) is 
along the real axis and it is bound between the chemical 
potentials /il and /ir. This is carried out over a numerical 
grid of sufficient fineness (dE)n- The calculation of p cq 
involves an unbound integral. This is performed over a 
coarse grid in the complex plane using a contour integral 
method^, since Gm is analytical^ and smooth in the 
imaginary energy plane. As shown in figure [SJa a number 
of numerical parameters must be chosen. First the lower 
limit of integration Eb must lie below the lowest lying 
molecular states and below the lowest electrode bands. 
Secondly, the poles of the Fermi functions (Matsubara 
frequencies) which lie within the contour must be taken 
into account. The integration is then performed by using 
Gaussian quadrature^ . 

The second self-consistent loop begins by calculating 
{G^, G^' } using equation (fl"4"|) . where the converged E H 
from the first loop is used and E F = 0. We then proceed 



and 






1.6 





1.6 


- 


3.2 





-3.2 








- 


1.6 



The c-p coupling a remains a free parameter. 

The parameters for the Au RC's match closely those 
used by Frederiksen o&^. The atoms in the leads are 
chosen as identical to the atoms in the scattering region, 
thus that a single onsite energy and hopping parameter 
characterise the electronic Hamiltonian. These parame- 
ters and the equilibrium potential p eq are chosen so that 
the differential conductance for the unperturbed system 
is Go = 2e 2 /h (perfect transmission). 



IV. RESULTS: H 2 -Pt JUNCTIONS 

A. self-consistent simulations 

The self-consistent G(V) in the range 0-200 meV for 
a ranging between and 3.1 eV/A are presented in fig- 
ure [SJ For this system the characteristic signature of e-p 
interaction is a drop in the conductance at a thrcshhold 
voltage Vthr- This signals the onset of inelastic electron 



7 



System parameter 


H 2 - Pt Au RC 




Symbol 


Value 


Value 


Units 


E F 


0.00 


0.00 


eV 


cm (molecule) 


-6.0 


0.0 


eV 


el (leads) 


0.00 


-1.00 


eV 


7l (leads) 


5 


-1.00 


eV 


7m (molecule) 


6.0 


-1.00 


eV 


7LM 


3.2 


-1.00 


eV 


T 


4.0 


4.0 


K 


m (atomic mass) 


1 


197 


amu 


K c 


21.82 


2.00 


eV/A 2 


K c _ c 


0.91 


1.00 


eV/A 2 


Numerical parameter 


Value 


Value 


Units 


[Vi»i,Vfinai] (Bias Range) 


[0,200] 


[0,31] 


meV 


Number of Bias points 


250 


120 


- 


(dE) F (Fock) 


0.1 


0.1 


meV 


Np (Fock) 


12600 


12600 


- 


rpi, (Hartree) 


4000 


4000 


- 


(dE) H (Hartree) 


.1061 


.0215 


meV 


cpi (Hartree) 


400 


200 


- 


cp2 (Hartree) 


400 


200 


- 


Eb (Hartree) 


-28.0 


-5.0 


eV 


Poles (Hartree) 


80 


80 




t-p (Tolerance) 


9.10" 8 


9.10" 8 


eV" 1 


in (Tolerance) 


1.10" 8 


1.10" 6 





TABLE I: Parameters used to simulate the H2-Pt junctions 
and Au RC's. The parameter K c is to the spring constant 
between the atoms in the chain, while K c _ e is to the spring 
constant between the molecule and the electrodes. The Fock 
grid and real Hartree grid are symmetric about Ef. 



processes involving the emission of phonons with ener- 
gies f2 w eVthr- We quantify this effect by denning the 
conductance drop (in units of Go) 

A thr = G(0) - G(V thr ) . (47) 




u.uj u.u ; j \j. i u u.v/j 

V [eV] V [eV] 



FIG. 6: Differential conductance and conductance drop at 
threshold for the H2-Pt junctions - V t h r is taken at 68.5 meV. 
Results obtained using the full SCBA are given for the dif- 
ferent initial conditions (panels a and c). In panels b) and 
d) we also show results obtained by setting E H = 0. In this 
second case Athr and G(V) agree well with the EST— (EST) 
up to a w 4.0 eV/A. The full SCBA disagrees with the EST 
at a considerably lower a « 1.8 eV/A. Above such a coupling 
strength G(V) depends on the initial conditions. 

fore the Fock one, the dependence on the initial condi- 
tions suggests that S H may be responsible for the be- 
haviour observed for a > a cr it- Stronger evidence to 
support this hypothesis is provided in figures [Hb and [SJi 
where the results for our second set of simulations in 
which S H is set to zero, are presented ("partial SCBA"). 
Figure [Hb shows good agreement between the partial 
SCBA and the EST to the much higher coupling of 
a ~ 4.0 eV/A. Moreover there is no evidence of any a cr it- 
This is confirmed by the G(V) curve obtained for a = 2.0 
eV/A and presented in figure [5Ji. 



Numerical simulations arc carried out in two differ- 
ent ways. First, we follow the exact numerical pro- 
cedure of figure^] ("full SCBA"), but we run simula- 
tions starting with different initial conditions, namely 
{r-Ajmit = + 0.003125 and {r A } init = - 0.003125. In 
figure Athr is plotted versus a demonstrating good 
agreement between the full SCBA and EST for low 
a. However, the two methods disagree for a beyond 
Q^crit ~ 1-8 eV/A. Athr peaks sharply above a cr it, be- 
yond which it becomes dependent on the initial condi- 
tion {r,\}init- This last situation is shown in figure Hjfc for 
a = 2.0 eV/A where two different G(V) curves are pre- 
dicted for different {rA}init and a low-bias conductance 
of 0.0125 Go is observed in stark disagreement with the 
unperturbed value of ~ 0.97 Gq. 

Since the Hartree self-consistent loop is performed be- 



B. Contribution from the individual modes 

We now analyze the origin the breakdown of the full 
SCBA. For V <C Vthr the inelastic current Ji ne i is strongly 
suppressed by Pauli exclusion principle. At low temper- 
ature (T=4.0 K) we can approximate the Fermi distribu- 
tions in J e i by step functions to obtain 

j/3(V) ~ tl I t(cj)6u. 
"-Jul 

We can now easily probe the contribution of the Hartree 
term to the conductance by considering the test Green 
function 

G$n(w) = [lu'Im - H u - S L - S R - rxR*}- 1 . (48) 
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This is used to evaluate the transmission coefficient at 
V = and it is useful to understand the influence of the 
individual modes A over the transmission. In figure [7] we 
show T(u>) as a function of r\ for the two longitudinal 
modes available in the H2-Pt system (see figure [8]). In 
the case of the rigid translational mode (A = 1 and fig- 
ure Ek) T(u>) is reduced in the region around E-p as |ri| 
increases. However the general shape of T(ui) is little af- 
fected. This is somehow expected from the shape of the 
matrix R 1 (eq. (05])), which indicates that mode 1 causes 
simply a change in the hopping parameters 7J2 and 734 
connecting the molecule to the leads. Importantly when 
is increased, 734 is reduced by the same amount and 
vice- versa. Moreover there is a symmetry T\ — ► — r\. 




co-E p [eV] 



FIG. 7: Transmission coefficients calculated using the Green 
function of Eq. (|48j for A = 1 (a) and A = 2 (b). The 
matrices R 1 and R 2 are fixed and a range of values for r\ has 
been chosen. 

The results for the symmetric mode A = 2 are pre- 
sented in figure [TJd. This time the peak in transmission 
is shifted in energy either to the left or to the right de- 
pending on the sign of r 2 . When shifted to the left, the 
peak is broadened while a shift to the right narrows it. 
In either case the transmission around Ep is reduced. 




FIG. 8: The longitudinal vibration modes of HVPt junctions. 
The A = 1 mode of energy 63 meV is the rigid translational of 
the H2 centre of mass, while the A = 2 mode is the symmetric 
mode of energy 432 meV. The atoms of the leads are fixed. 

From the figures [6J:, [7^,, and[7}D one can conclude that 
the deviation of the zero bias differential conductance 
from its unperturbed (no e-p interaction) value is a mea- 
sure of the magnitude of the e-p perturbation. We define 



this deviation as 

AG = Go - G(0) = G - T(E F ) , (49) 

where the last equality is valid for T — > 0. Figures [5^ and 
12b show the estimated deviation (in units of Go) versus 
the weighting coefficients r\. The maximum deviation, 
AG = I, occurs when the chain actually breaks as for 
mode I and r\ = 1. As expected, the curve in figure^ 
for mode 1 is symmetric about while figure [SJa for mode 
2 is not. 




[arb. units] 



FIG. 9: Estimate of the deviation AG caused by an individual 
mode A as a function of r\. The insets shows the region of 
small perturbation. 



C. Discussion of the self-consistent results 

We finally re-analyze our self-consistent results in the 
light of the discussion in the previous section. Figure 
[TOh andflOb show the self-consistently calculated r\ as a 
function of a for V = 0, while figure [TUb shows AG as de- 
fined in equation (T4"9")) . The critical point a cr it ~1.8 cV/A 
marks a sharp transition in the behaviour of the weight- 
ing coefficients. This is evident in the abrupt change 
of magnitude and behaviour of AG (in figure [TUb) . In 
fact while for a < a cr n the full SCBA agrees well with 
the EST, the two differ sharply as soon as a > a cr it- 
Going into more details we note that r% is identically 
zero (~ 10~ 14 ) for a < a CT it, while r 2 is small and neg- 
ative. Importantly both n and ri are independent of 
{r,\}imt- By contrast, for a > a cr it, r-i remains indepen- 
dent of {r.\}imt but this is not the case for r\. In fact 
we obtain a strong dependence over the initial condi- 
tions with positive (negative) T\ for positive (negative) 
{?~A}mit- Importantly none of these features are found in 
the case when we neglect the Hartree self-energy. 

Since |ri| is two orders of magnitude larger than |r*2 1 , 
it will largely determine E H . As figure [TOh shows, 7*1 
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FIG. 10: Converged self-consistent values of r\ and deviation 
AG (Eq. (|31|) as a function of the e-p coupling strength a for 
V = 0. A transition is apparent at a = 1.865 ± 0.005 eV/A 
for the full SCBA. Panels a) and b) show respectively n and 
T2 for the full SCBA for different initial conditions, while in 
c) results for the partial SCBA and EST are also included. 



varies roughly linearly with a above a cr it, therefore the 
self-consistent AG follows the estimated curve of figure 

in this region of a. The magnitude of n for a > a CT n 
suggests that the interaction with phonons becomes a 
strong perturbation of the electronic system. We define 
the region a < a cr it as the weak coupling regime and the 
region a > a cr it as the strong coupling regime. 

Figure [TO] adequately explains the causes of the mas- 
sive reduction of G(V) with respect to its unperturbed 
value observed in figured]; at V = 0. Notably as figure 
[f^i shows, A t hr calculated with the full SCBA starts de- 
viating from the EST result for a ~ 1.75 eV/A, i.e. at a 
value lower that a cr it = 1.865±0.005 eV/A calculated for 
V = 0. This seems to suggest that the critical value of 
a for the breakdown of the full SCBA somehow depends 
on the bias. Moreover at finite bias a cr it is character- 
ized by a peak in Athr (a) for positive {r\\}; n it and by a 
discontinuity for negative {r,\}init (see figure^). 

In order to explore the onset of the breakdown of the 
SCBA at finite bias in figure [TTk we present G(V) for 
a = 1.84 eV/A i.e. just below the zero-bias critical value. 
In addition, in figure [Tib we plot the dominant coefficient 
r\ for a small range of a about a cr it at three different 
bias, V = 0, V = 0.1 and V = 0.2 Volt. A clear result 
from figure [TTb is that the presence of S H introduces 
a reduction of G(V) with bias not present at V = 0. 
For a < a cr it figure [Tib shows that \r%\ is a function of 
bias whose value increases as the bias increases. The 
difference between r% at V = and at finite V explains 
the deviation below cv cr it and also the origin of the peak 
above it. 

The discontinuity in Athr for {r,\}i n it > of figure^ is 
explained by the discontinuities observed in r\ (Fig. II lb ") 
and in the current J (V) (Fig. [TTbV In fact at V = 0, 



a) 1 

0.975 
0.95 
:= 0.925 

0.875 
0.85 
0.825 

c) 




a = 1.84 eV/A 



4— < V=0.1eV 
►— ► V =0.2 eV 

V = (1^,= +0.003125 
»-oV = 0{i,L= -0.003125 








0.1 
V [eV] 
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FIG. 11: Behaviour of the SCBA in the region above and be- 
low the zero-bias « C rit- For a < a cr it, G(V) for the full and 
partial SCBA are compared in panel a), showing that the con- 
tribution arising from E H is bias-dependent. For a > a cr it, 
J 13 calculated with the full SCBA is presented in c) for two 
different initial conditions. In panel b) the functional depen- 
dence of the weighting coefficient n upon bias is investigated 
for a range of coupling strengths and bias. Finally, in panel 
d) we show the magnitude |ri| obtained from two full SCBA 
simulations with different initial conditions and a > a cr it- 



ri is single valued as long as a < a cr it- For a > a cr it 
instead r\(a) has a parabolic shape symmetric about 
r\ = 0: the {r^linit < solution traces out the lower arm 
of the parabola and the {rAjinit > solution follows the 
upper arm as a increases. For V > 0, it is seen that 
the two solutions for r\ are identical, and asymptotically 
approach the lower arm of the V = curve from be- 
low for a > ev cr it. However, as a is further increased the 
{''Ajinit > solution jumps discontinuously above zero 
and then asymptotically approaches the upper arm, again 
from below. 

The discontinuity in n(a,V) is determined by both 
the bias and the initial conditions. Generally, it is 
found that such discontinuity occurs for lower bias first; 
r\ jumps discontinuously for V = 0.1 Volt before it 
does for V = 0.2 Volt. This explains the behaviour of 
the {r A } ini t = +.003125 solution for J?(V) in figured]]; 
which leads to a peak in its derivative (the conductance 
G(V)) and explains the discontinuity of A t hr- 

We note that after the discontinuity in n for V > 0, 
the two solutions are no longer symmetric about n = 
and do not converge to the V = solutions of cither 
arm until a « 3.0 eV/A. This is highlighted in figure [TTB 
where for the two solutions is plotted in a range of 
a just above a cr ;t at V = 0.2 Volt. Such differences 
explain why A t hr is not independent of the initial con- 
ditions beyond the discontinuity and also the different 
curves observed for G(V) in figure [HJi. 
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Mode 1(R) 



Mode 3 (A) 



• • II 



a [eV/Al 



FIG. 12: G(V) for the 4C illustrating the onset of in- 
elastic processes at threshhold voltages. These are associ- 
ated to the symmetric modes of energy Q2 = 6.5 meV and 
fU = 12.2 meV. The rigid and anti-symmetric modes of en- 
ergies respectively =3.1 meV and Q3 = 9.8 meV have 
no effect for chains comprising an even number of atoms. No 
overall shift in G(V) is observed as a lies below the critical 
coupling a cr it which is clearly determined from the inset in 
panel b. 



We make a final comment about the discontinuities 
seen in figure [Tib . The value of a at which these occur is 
dependent on the initial conditions as mentioned already. 
Thus for a particular {r>}i n i t it may be possible to reach 
the upper solution at a cr it for all bias, so that r\ has a 
parabolic shape for bias while being asymmetric about 0. 
We have not observed this and regard a cr it as uniquely 
defined for V = only. 



V. RESULTS: Au Chains 

By using the procedure outlined in section ITTTI and the 
parameters of table |TJ the full SCBA is used to calculate 
the transport of the Au RC's. In general we observe 
a behaviour similiar to that of the E^-Pt system. For 
3C, 4C, 5C, and 6C a weak coupling regime is identified 
where the shift AG and the weighting coefficients are 
zero. A critical coupling strength a CI1 t for V = was 
discovered for each of the chains investigated with values 
a crit « 0.85,0.9,0.82,0.83 eV/A respectively for 3C, 4C, 
5C, and 6C. For weak coupling AG matches closely the 
values calculated in previous work a 24 i 35 . As an example, 
in figure [T2l we present G(V) and AG (a) for 4C. It is seen 
that the modes symmetric with respect to the centre of 
the scattering region (even numbered) induce drops in 
G(V) at threshhold voltages corresponding to the energy 
of the modes. In general, only the symmetric modes are 
active in RC's containing an even number of atoms and 
conversely the rigid and anti-symmetric modes arc active 



Mode 2(S) 

« • W • 



Mode4(S) 



9) ww 



FIG. 13: Vibration modes of Au 4C. Modes 2 and 4 are sym- 
metric modes (S) about the centre of the scattering region. 
Mode 1 is the rigid translational mode (R) while mode 3 is 
anti-symmetric (A). For all chains considered, the mass of the 
lead atoms are taken sufficiently larger than that of the chain 
atoms so that they do not vibrate. 
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FIG. 14: Test of the convergence of the weighting coefficient 
T4 for the 4C and a > a cr it using a range {rA}mit where n is 
the number of iterations. Results for the calculations with 
positive starting values are indicated with dashed lines while 
those with negative initial values are shown with solid lines. 
A single negative and positive solution for the converged 14 is 
found. 



for an odd RC. 

The transition from weak to strong coupling can be ap- 
preciated for the 4C by looking at the inset of figure [12b. 
where for a > a cr it two different initial conditions lead to 
two different AG. The 6C shows the same behaviour of 
the 4C, while the 3C and 5C show a single curve for AG 
due to the symmetry of the weighting coefficients in the 
strong regime. Beyond ev cr it, for all the RC's simulated 
G(V) is reduced to zero as a increased and the shift AG 
rises to Go . 

Finally we want to investigate further the existence of 
multiple solutions depending on the initial conditions in 
the self-consistent loops. As an example of how conver- 
gence is achieved in figure [14] we show the coefficient 
as a function of the iteration number n for the 4C plotted 
for a sing le bias V = 0.02 Volt and coupling 0.902 cV/A 
(a > a cr ;t). A number of simulations were run with dif- 
ferent initial conditions. The figure indicates the exis- 
tence of two stable minima: simulations which start at 
{fAj-init > (dashed lines) converge to the same positive 
final value, while simulations initialised with {rAjinit < 
converge to the same r$ < value. We note that the min- 
ima are not symmetric about rg = 0. The solution r\ = 
appears to be a minima in the weak coupling regime, but 
becomes unstable and evolves to a local maximum in the 
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FIG. 15: Convergence test. The zero-bias r& is plotted versus 
the iteration number n for {rA}i n it = —0.09 and a > a cr it- 
Here the lower bound for the integration of the charge density 
p, EB, is varied. The lower band-edge of the leads lies at - 
2 eV. Clearly when EB > —2 eV r§ does not converge fully, 
i.e. it depends on EB. 
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FIG. 16: T(u>) for a 6C simulation where only the sixth 
mode of energy — 12.6 meV is included. The curves 
show results obtained when E H is calculated using the Simp- 
son's rule ({rgjinit = in the legend) or the contour method 
({r 6 }init = -0.00975 in the legend). 



strong coupling regime. 

In figure \T5\ i'6 is plotted versus the number of itera- 
tions n for a > a cr it- This time we run different simula- 
tions in which the lower bound of the energy integration 
grid, EB, is changed. In particular we explore situations 
where EB is not below the lower band-edge of the leads, 
that for our choice of parameters lies at -2 eV (see table 
QJ . The figure indicates that if EB > — 1 eV r$ converges 
to zero so that the sixth mode gives no contribution to 
S H . For —2 < EB < —1 r§ is nonzero, however the con- 
verged value differs for EB= -1.5 eV and EB= -2.0 eV, 
i.e. it is sensitive on the grid lower bound. Finally for 
EB < — 2.0 eV the bands of the leads are entirely included 
in the integral and converges to a value of approxi- 
mately —0.25 eV which is independent of our choice of 
EB. From this simple analysis it appears that cutting the 
integration grid can result in the erroneous suppression 
of the Hartree self-energy, i.e. in a drastic underestimate 
of its contribution. This produces a fortuitous suppres- 
sion of the SCBA breakdown, since the agreement be- 
tween SCBA and EST is usually improved when S H is 
neglected. 

Finally we test the robustness of our integration 
method. Figure [16] shows the transmission coefficients 
T(uj) for a 6C at a bias of V = 1 mV where only the 
sixth mode is considered. We run two simulations. In 
the first case, the full SCBA is used with the numer- 
ical parameters taken from table |T] and initial condi- 
tion {rgjinit = —0.00975. In the second case, the in- 
tegration method outlined in section IIII Bl is replaced 
by Simpson's rule along the real axis. The integra- 
tion range used is [-3.093,0.1] eV with a grid fineness 
(dE) H = 2.129 x lO -5 eV. The Hartree self-consistent 
loop is started with {ro\i n it = 0.0 and finished when 
c" < tjj = 10~ 12 (this tolerance was also used in the first 
case). All other parameters of the calculations are iden- 
tical in the two cases. The figure shows that the two 
numerical methods produce the same T(uj), specifically 
in the range of applied bias shown in the inset. 



VI. CONCLUSIONS 

By using a simple ID tight binding model we have 
investigated the breakdown of the SCBA as a function 
of the e-p coupling strength a. We have identified two 
regimes. In the weak coupling regime there is a unique 
solution for both E H and S F , independently of the ini- 
tial conditions. In particular the Hartree self-energy is 
small and has little effect on the final conductance. In 
this weak coupling regime the characteristic conductance 
drops at voltages corresponding to the various phonon en- 
ergies compare well with those calculated with the EST 
method. 

As the coupling parameter a is increased beyond some 
critical value a cr it a sharp transition to the strong cou- 
pling regime occurs. In this limit the self-consistent E H 
becomes unstable with respect to the initial conditions 
and exhibits multiple values for the same voltage. This 
results in a conductance that sharply deviates from that 
obtained with the exact EST. Such a sharp transition 
suggests a breakdown of the SCBA, indicating that the 
electron-phonon interaction cannot be treated perturba- 
tively for a > a CT it- Interestingly such a breakdown 
is suppressed when the Hartree self-energy is neglected 
completely from the calculation, as routinely done in 
practice. Our results however set a warning to quantita- 
tive calculations based on e-p parameters extracted from 
density functional theory. For these, no information is 
available on whether or not the obtained c-p coupling 
strength is cither below or above the critical value for 
the SCBA to breakdown. Therefore, strictly speaking, 
one rarely knows whether the SCBA is applicable. 
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